The manipulation of light propagation by means of periodic photonic structures is an important scientific idea, which has great values for potential applications in various branches of engineering. In particular, arrays of evanescently coupled waveguides are unique structures that exhibit the peculiar properties. The linear and nonlinear properties of uniform waveguide arrays have been studied theoretically and experimentally, such as anomalous discrete diffraction, Bloch oscillations, and discrete solitons.^[1] In the last few years, the investigation of nonuniform waveguide arrays has received increasing attention, since a more complex engineering of the periodic structure can provide further degrees of freedom. For example, the zigzag waveguide arrays were proposed to obtain diffraction management.^[2] Moreover, binary arrays composed of waveguides with alternating widths and spacings have been thoroughly studied to exhibit double refraction and gap solitons.^[3–5] In the sub-wavelength regime, the plasmonic waveguides provide a new way to break through the diffraction limit. However, plasmons in metals exhibit relatively strong losses and cannot easily be controlled. Recently, the graphene-supported surface plasmon polaritons (SPPs) have drawn a great deal of attention due to their unique physical properties.^[6,7] Graphene exhibits a stronger SPP field confinement and lower propagation loss compared with metals.^[8] Moreover, the surface conductivity of graphene can be flexibly tuned by chemical doping or external static electric and magnetic fields. The above estimation highlights the potential of graphene plasmonics in new tunable optical devices at nanoscale.^[9,10]

In this work, we shall comprehensively study the diffraction properties of the graphene sheet arrays (GSAs) composed of two types of graphene sheets with different chemical potentials. The theoretical calculations are based on the transfer matrix method (TMM)^[11–13] and coupled-mode theory (CMT). Furthermore, numerical simulations are performed to assess the validity of the analytical treatment.

The system consists of N-layer graphene sheets, which is shown in Fig. 1. The position of graphene sheets is denoted by x _n , where , 0, d, …. The materials separated by graphene sheets are dielectric and have the permittivity of . The graphene surface conductivity is modeled by the Kubo formula. When the chemical potential of graphene is much larger than the interacting photon energy , the surface conductivity of graphene has a Drude-like form ,^[14] where e is the electron charge and τ is the momentum relaxation time of electrons. The graphene sheets possess different chemical potentials labeled as and ( ), while the chemical potential difference between sheets A and B is . The relaxation time is set as 1 ps at room temperature T = 300 K.^[15] In our study, the incident wavelength and the interlayer space d = 70 nm are initially considered. Here the interlayer space is large enough, which suggests the weak coupling of adjacent graphene sheets in the array.^[16]

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Schematic of binary GSAs, consisting of two interleaved graphene waveguides A and B with different chemical potential and , equally spaced by d.

We only consider TM polarized SPPs propagating along the z direction, and the magnetic field in the region between the position of and d can be written as the superposition of waves in opposite directions:

(1)

According to Maxwell's equations, the tangential electric field reads as

(2)

(3)

(4)

To analyze the existence conditions of the surface modes, we describe their formations based on coupled-mode theory. By considering only nearest-neighbor couplings, the normalized amplitudes of field distributions in binary GSAs obey the coupled-mode equations^[17,18]

(5)

(6)

The diffraction curves in the binary GSAs for d = 70 nm, as eV are plotted in Fig. 2. The binary GSAs support two minibands separated by a bandgap . In contrast with the metal-dielectric arrays, the bandgap is determined at the edge of the Brillouin zone ( ). The curves obtained from Eq. (6) (blue curves) agree with the TMM method (red curves), which verifies that the coupled-mode model holds well in the binary GSAs structure.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Diffraction relations of the binary GSAs for d = 70 nm and as , comprising two bands separated by 2δ. The red and blue curves are the results calculated by the TMM method (Eq. (4)) and coupled-mode equation (CME, Eq. (6)), respectively.

The diffraction curves in the binary GSAs for different as are plotted in Fig. 3. For , there is no gap between two symmetric minibands (as shown by the dotted curves). While for a finite difference of chemical potentials, e.g., , a gap opens at the edge of the Brillouin zone. As shown in Fig. 3, when is changed from 0 to 0.012 eV, the bandgap enlarges as increases. In the case of discrete excitation with nonzero , there should be two propagating SPP fields inside the array with different diffraction relations .

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Diffraction relations of the binary GSAs in the reduced Brillouin zone for d = 70 nm with different initial chemical potential difference between graphene sheets A and B ( , , and eV), as is fixed at 0.15 eV.

The steady electromagnetic field distributions are numerically calculated to verify the above theoretical predictions. The computations are performed by using the finite-difference frequency-domain (FDFD) method,^[19,20] where the graphene sheets are assumed to be freestanding in air with . Graphene is equivalent to a very thin metal film with a thickness of . Thus, we can define an equivalent bulk permittivity for graphene given by while the air impedance . The minimum mesh size equals 0.2 nm in the FDFD calculation. To achieve the two-way splitting of the input power in the binary GSAs effectively, a TM-polarized Gaussian beam with a small tilt is incident into the binary GSAs. Under this condition, the incident Bloch momentum is 0.25π. The full width at half-maximum (FWHM) of the incident Gaussian beam is about 0.21 m, which excites mostly three waveguides. The distributions of simulated electric field intensity ( ) of SPP beams at different are shown in Figs. 4(a) and 4(b), and the input and output field intensity distribution is shown in Fig. 4(c). When , the initial tilted input beam experiences discrete diffraction in the homogeneous GSAs. When , the beam splits into two paths due to different refraction angles of the two diffraction bands. The refraction angle could be expressed as . According to the diffraction curves, the refraction angles have two different values for . In this case, the normalized transverse intensity distributions at two different propagation distances, and m, are also depicted in Fig. 4(c), respectively. Here the initial SPP beam ( m) splits into two beams belonging to the two bands of the array at the propagation distance of m. The simulated results of the SPP propagation behaviour coincide with our theoretical predictions.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) Steady electric field ( ) distributions of SPPs (a) in the homogeneous GSAs with and (b) in the binary GSAs with , . (c) The input (blue line) and output (red line) field intensity from the FDFD simulation.

The diffraction properties of the binary GSAs are investigated by theoretical analyses and performing FDFD simulations. The pair of graphene sheets with different chemical potentials are stacked alternatively. Two diffraction bands are separated by a bandgap at the edge of the Brillouin zone in these binary GSAs, which are very sensitive to the chemical potential difference between adjacent graphene sheets. The splitting of the SPP beams in these GSAs has been found, which is due to the two different diffraction bands. These surface modes in binary GSAs could be used to manipulate SPPs in the nano-optics devices.